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The direction-of-arrival (DOA) estimation problem involves the localization of a few sources from a 
limited number of observations on an array of sensors, thus it can be formulated as a sparse signal re¬ 
construction problem and solved efficiently with compressive sensing (CS) to achieve high-resolution 
imaging. On a discrete angular grid, the CS reconstruction degrades due to basis mismatch when the 
DOAs do not coincide with the angular directions on the grid. To overcome this limitation, a con¬ 
tinuous formulation of the DOA problem is employed and an optimization procedure is introduced, 
which promotes sparsity on a continuous optimization variable. The DOA estimation problem with 
infinitely many unknowns, i.e., source locations and amplitudes, is solved over a few optimization 
variables with semidefinite programming. The grid-free CS reconstruction provides high-resolution 
imaging even with non-uniform arrays, single-snapshot data and under noisy conditions as demon¬ 
strated on experimental towed array data. 


PACS numbers: 43.60.Pt, 43.60.Jn, 43.60.Fg 


I. INTRODUCTION 

Sound source localization with sensor arrays involves 
the estimation of the direction-of-arrival (DOA) of (usu¬ 
ally a few) sources from a limited number of observations. 
Compressive sensin^i^^ (CS) is a method for solving such 
underdetermined problems with a convex optimization 
procedure which promotes sparse solutions. 

Solving the DOA estimation as a sparse signal re¬ 
construction problem with CS, results in robust, high- 
resolution acoustic imaginJ^HSl^ outperforming traditional 
methodd^ for DOA estimation. Furthermore, in ocean 
acoustics, CS is shown to improve the performance of 
matched field processin^^, which is a generalized beam¬ 
forming method for localizing sources in complex envi¬ 
ronments (e.g., shallow water), and of coherent passive 
fathometry in inferring the number and depth of sedi¬ 
ment layer interfacesISl. 

One of the limitations of CS in DOA estimation is basis 
mismatcff^ which occurs when the sources do not coin¬ 
cide with the look directions due to inadequate discretiza¬ 
tion of the angular spectrum. Under basis mismatch, 
spectral leakage leads to inaccurate reconstruction, i.e., 
estimated DOAs dev iating from the actual ones. Em¬ 
ploying finer gridd^^ alleviates basis mismatch at the ex¬ 
pense of increased computational complexity, especially 
in large two-dimensional or three-dimensional problems 
as encountered in seismic imaging for examplel^"!^. 

To overcome basis mismatch, we formulate the DOA 
estimation problem in a continuous angular spectrum 
and introduce a sparsity promoting measure for general 
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signals, the atomic norirfi^. The atomic norm mini¬ 
mization problem, which has infinitely many unknows, 
is solved efficiently over few optimization variables in the 
dual domain with semidefinite programmin^i^. Utiliz¬ 
ing the dual optimal variables, we show that the DOAs 
are accurately reconstructed through polynomial root¬ 
ing. It is demonstrated that grid-free CS gives robust, 
high-resolution reconstruction also with non-uniform ar¬ 
rays and noisy measurements, exhibiting great flexibility 
in practical applications. 

Polynomial rooting is employed in several DOA estima¬ 
tion methods to improve the resolution. However, these 
methods involve the estimation of the cross-spectral ma¬ 
trix hence they require many snapshots and stationary 
incoherent sources and are suitable only for uniform lin¬ 
ear arrays (ULA)P^. Grid-free CS is demonstrated not to 
have these limitations. 

Finally, we process acoustic datcP^from measurements 
in the North-East (NE) Pacific with grid-free CS and 
demonstrate that the method provides high-resolution 
acoustic imaging even with single-snapshot data. 

In this paper, vectors are represented by bold lower¬ 
case letters and matrices by bold uppercase letters. The 
symbols ^ denote the transpose and the Hermitian 
(i.e., conjugate transpose) operator respectively on vec¬ 
tors and matrices. The symbol * denotes simple conju¬ 
gation. The generalized inequality X ^ 0 denotes that 
the matrix X is positive semidefinite. The £p-norm of a 

vector X e C” is defined as ||x||p = 
extension, the .^o-nonn is defined as ||x||o = 

The paper makes heavy use of convex optimization the¬ 
ory; for a summary see App. 
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II. DISCRETE DOA ESTIMATION 

The DOA estimation problem involves the localization 
of usually a few sources from measurements on an array 
of sensors. For simplicity, we assume that the sources 
are in the far-held of the array, such that the wave- 
held impinging on the array consists of a superposition 
of plane waves, that the processing is narrowband and 
the sound speed is known. Moreover, we consider the 
one-dimensional problem with a uniform linear array of 
sensors and the sources residing in the plane of the array. 

The location of a source is characterized by the di¬ 
rection of arrival of the associated plane wave, 9 S 
[—90°,90°], with respect to the array axis. The propa¬ 
gation delay from the ith potential source to each of the 
M array sensors is described by the steering (or replica) 
vector, 

a(6l,) (1) 

where A is the wavelength and d is the intersensor spac¬ 
ing. 

Discretizing the half-space of interest, 9 G [—90°,90°], 
into N angular directions the DOA estimation problem 
is expressed in a matrix-vector formulation, 

y = Ax, (2) 

where y G is the vector of the wavefield measure¬ 
ments at the M sensors, x G C'^ is the unknown vector 
of the complex source amplitudes at all N directions on 
the angular grid of interest and A is the sensing matrix 
which maps the signal to the observations, 

AmxAt = [a(0i), • • • , a(0Ar)]. (3) 

In the presence of additive noise n G C^, the measure¬ 
ment vector is described by, 

y = Ax + n. (4) 

The noise is generated as independent and identically 
distributed (iid) complex Gaussian. The array signal-to- 
noise ratio (SNR) for a single-snapshot is used in the sim¬ 
ulations, defined as SNR=20log^o (|jAx|j 2 /||n|j 2 ), which 
determines the noise .^ 2 -norm, ||n ||2 = |j Ax|| 2 l 0 “®^^/^°. 

A. Sparse signal reconstruction 

Practically, we are interested in a fine resolution on the 
angular grid such that M < N and the problem ([^ is 
underdetermined. A way to solve this ill-posed problem is 
to constrain the possible solutions with prior information. 

Traditional methods solve the underdetermined prob¬ 
lem ([^ by seeking the solution with the minimum £ 2 - 
norm which fits the data as described by the minimiza¬ 
tion problem, 

min ||x|j 2 subject to y = Ax. (5) 

xGC« 

The minimization problem ([^ is convex with analytic 
solution, X = (AA^) ^ y. However, it aims to min¬ 

imize the energy of the signal rather than its sparsity, 
hence the resulting solution is non-sparse. 


Conventional beamformin^^ (CBF) is the simplest 
source localization method and it is based on the £ 2 - 
norm method with the simplifying condition AA^ = 1m- 
CBF combines the sensor outputs coherently to enhance 
the signal at a specific look direction from the ubiquitous 
noise yielding the solution, 

XCBF = A^y. (6) 

CBF is robust to noise but suffers from low resolution 
and the presence of sidelobes. 

A sparse solution x is preferred by minimizing the £ 0 - 
norm leading to the minimization problem, 

min ||x|jo subject to y = Ax. (7) 

xGC” 

However, the minimization problem 0 is a non-convex 
combinatorial problem which becomes computationally 
intractable even for moderate dim ensions. The break¬ 
through of compressive sensin^i^^ (CS) came with the 
proof that for sufficiently sparse signals, K « N, 
K < M, and sensing matrices with sufficiently incoher¬ 
ent columns the minimization problem Q is equivalent 
to the minimization problem, 

min ||x||i subject to y = Ax, (8) 

xGC” 

where the £o“norm is replaced with the £i-norm. The 
problem (|^ is the closest convex optimization problem 
to the problem Q and can be solved efficiently by convex 
optimization even for large dimension^^. 

For noisy measurements Q, the constraint in Q be¬ 
comes ||y — Ax ||2 < e, where e is the noise floor, i.e., 
Iln ||2 < £■ Then, the solution i^, 

xcs = argmin||x||i subject to ||y — Ax II 2 < e, (9) 

xGC« 

which has the minimum fi-norm while it fits the data up 
to the noise level. 

Herein, we use the cvx toolbox for disciplined convex 
optimization which is available in the Matlab environ¬ 
ment. It uses interior point solvers to obtain the global 
solution of a well-defined optimization problercPHl. Inte¬ 
rior point methods solve an optimization problem with 
linear equality and inequality constraints by transform¬ 
ing it to a sequence of simpler linear equality constrained 
problems which are solved iteratively with the Newton’s 
method (iterative gradient descent method) increasing 
the accuracy of approximation at each stepP^. 

B. Basis Mismatch 

CS offers improved resolution due to the sparsity con¬ 
straint and it can be solved efficiently with convex opti¬ 
mization. However, CS performance in DOA estimation 
is limited by the coherence of the sensing matrix A (see 
Refill, described by t he re stricted isometry propertjP^, 
and by basis mismatc H^^F^ I pg inadequate discretiza¬ 
tion of the angular grid. Herein, we demonstrate a way 
to overcome the limitation of basis mismatch by solving 
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FIG. 1. (Color online) CS performance in DOA estimation in 
terms of the discretization of the angular space. A standard 
ULA is used with M=8 sensors, d/X =1/2 and SNR = 20 dB. 
CBF and CS (*) reconstruction of two sources (o) (a) at 0° 
and 15° on a grid [—90°:5°:90°], (b) at 0° and 17° on a grid 
[—90°:5°:90°] and (c) at 0° and 17° on a grid [—90°:1°:90°]. 


and the measurement vector of the sensor array is, 


yMxl=^MX, ( 12 ) 

where J^m is a linear operator (inverse Fourier transform) 
which maps the continuous signal x to the observations 
y e C^. 

In the presence of additive noise, n S , the mea¬ 
surement vector is described by, 

y = J^MX -b n, (13) 


similarly to Eq. Q. 


the £i-minimization problem on a grid-free, continuous 
spatial domain. 

The fundamental assumption in CS is the sparsity of 
the underlying signal in the basis of representation, i.e., 
the sensing matrix A. However, when the sources do not 
match with the selected angular grid, the signal might 
not appear sparse in the selected DFT basiJ^il. Figure 
shows the degradation of CS performance under basis 
mismatch due to inadequate discretization of the DOA 
domain in FFT beamforming. 

To increase the precision of the CS reconstruction, 
Malioutov et alW and Duarte and Baraniul^i^ propose an 
adaptive grid refinement. The adaptive grid refinement 
aims at improving the resolution of CS reconstruction 
without significant increase in the computational com¬ 
plexity by first detecting the regions where sources are 
present on a coarse grid and then refining the grid locally 
only at these regions. Grid refinement is an intuitive way 
of circumventing basis mismatch. However, the problem 
of basis mismatch is avoided only if the problem is solved 
in a continuous setting, particularly for moving sources. 


IV. GRID-FREE SPARSE RECONSTRUCTION 


To solve the underdetermined problem (12) (or equiv¬ 
alently the problem (HI) in favor of sparse solutions, 
we describe an optimization procedure which promotes 
sparsity on a continuous optimization variable. 


A. Atomic norm 


In the discrete formulation (|^ of the DOA estimation 
problem, the prior information about the sparse distri¬ 
bution of sources is imposed through the £i-norm of the 
vector X to obtain sparse estimates ([^ . By extension, in 
the continuous formulation (131, we introduce the atomic 
nornP^, ll’H^, as a sparsity promoting measure for the 


continuous signal x{t) in Eq. (10) defined as. 


K 


|2:|U 


= Ei' 


(14) 


III. CONTINUOUS DOA ESTIMATION 

In the continuous approach, the A-sparse signal, x, is 
expressed as. 


K 

x{t) ='^XiS{t - ti), (10) 

i=l 

where S C is the complex amplitude of the 1th source, 
ti = sin 9i is its support, i.e, the corresponding DOA, on 
the continuous sine spectrum T = [—1,1] (with T C T 
the set of the DOAs of all K sources) and 6{t) is the 
Dirac delta function. 

The sound pressure received at the mth sensor is ex¬ 
pressed as a superposition of plane waves from all possible 
directions on the continuous sine spectrum T, 


r ^ 


In other words, the atomic norm is a measure for contin¬ 
uous signals equivalent to the ^i-norm (which is defined 
only on vector spaces). Hence, the atomic norm is a con¬ 
vex function which promotes sparsity in a general frame¬ 
work. Eor a discrete grid the atomic norm corresponds 
to the £i-norm. 

To clarify the analogy between the .^i-norm and the 
atomic norm and justify the term atomic, consider that 
the vector x G can be interpreted as a linear combi¬ 
nation of N unit vectors. The unit vectors, in this case, 
are the smallest units, or atoms, in which the vector x 
can be decomposed into. The £i-norm is the sum of the 
absolute values of the weights of this linear combination 
of atomd^. 

Analogously, the continuous signal (10) can be in¬ 
terpreted as a linear combination of K delta functions 
S{t — ti), serving as atoms for the continuous signal x{t) 
and the atomic norm is the sum of the absolute values of 
the weights of the linear combination of these atomJ^^. 
Even though there are infinitely many atoms in the con¬ 
tinuous case, only few of those, K < M, constitute the 
signal and the sum in (14) is finite. 
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B. Primal problem 

Utilizing the convex measure of the atomic norm, the 
DOA estimation in the continuous angular space is solved 
with the sparsity promoting minimization problem, 


min||a;||^ subject to y = J-mX. 


(15) 


Since the optimization variable x is a continuous pa¬ 


rameter, the primal problem (15) is infinite dimensional 


and cannot be solved as such. It is possible to approx¬ 
imate the continuous variable a: on a discrete grid and 
solve the £i-norm optimization problem Q. This would 
increase the computational complexity significantly when 
the discretization step is reduced to improve precision. 
An alternative to this, is to gradually refine the dis¬ 
cretization stepP. However, we show that by solving the 
dual problem instead, there is no need to employ a dis¬ 
crete approximation of the continuous variable, x. 


C. Dual problem 


and cos (j)i < 1. Thus, for Xi ^ 0, is a unit vector 

in the direction of x,-, 


i — Xij\xi\^ Xi ^ 0 


< 1 , 


Xi = 0 . 


( 20 ) 


Maximizing the dual function (191 constitutes the dual 
problem, 

max Re Tc^yl subject to ||J^^c||oo < 1- (21) 

cgC" 


Since the primal problem (15) is convex with linear equal¬ 


ity constraints (All), strong duality holds assuring that 


the maximum of the dual problem (21) is equal to the 
minimum of the primal problem. 


The dual problem (21) selects a vector c G which is 
maximally aligned with the measurement vector y G 
while its beamformed amplitude is bounded by 

unity across the whole angular spectrum. At the angu¬ 
lar direction corresponding to the DOA of an existing 
source, the beamformed dual vector (201 is equal to the 
normalized source amplitude. 


To formulate the dual problem to the problem (15) (see 
Appendix for details), we construct the Lagrangian 
by making the explicit equality constraints, y = ^mx, 
implicit in the objective function, 


L{x, c) = ||x|U + Re [c^ (y - J^mx)] 


(16) 


where c G is the vector of dual variables. 

The dual function g[c) is the infimum, i.e., the greatest 
lower bound, of the Lagrangian, L(x,c), over the primal 
optimization variable x. 


D. Dual problem using semidefinite programming 


The dual problem (21) is a semi-infinite programming 


problem with a finite number of optimization variables, 
c G and infinitely many inequality constraints, 

which is still intractable. 

Define the dual polynomial, 


M-l 


H{z) = J-^c = ^ 


CmZ = 


m—Q 


M-l 

E 

771=0 


( 22 ) 


g{c) = inf L(x, c) 

X 

= Re [c^y] -f inf (||x|U - Re [c^J^mx]) . 


(17) 


To evaluate the second term in ( |l7| ) we note that 

H 

for every x*. Re [(c^ J'm) . x^] = Re 
|(j^j(|c) .||xi| cos(j)i, where 4>i is the angle between Xi and 
Then, 


|xJ-Re 




mC 


. H 


= \x^\[l-\{j^^c)^\cos (jji] 
> |x,|[l-|(J'^c)J] . 


(18) 


Note that is a trigonometric polynomial (B11, of the 
variable z{t) = t G T, with the dual variables 

c = [co, • • • , cm-i]^ as coefficients and degree M — 1. 
The inequality constraint in Eq. (21) implies that the 


dual polynomial has amplitude uniformly bounded for all 
t G T; see Eq. (B7|. Making use of the approximation 
in Eq. (B6) for bounded trigonometric polynomials, the 


constraint in Eq. (211 can be replaced with finite dimen¬ 
sional linear matrix inequalities. Thus, the dual prob¬ 
lem is solved with semidefinite programmin^^siSi]^ j ^ 
convex optimization problem where the inequality con¬ 
straints are linear matrix inequalities with semidefinite 
matrices, 


The lower bound in ( [I^ is nonnegative if |.F^c| is less 
than one, max^j J < 1, and the infimum is zero. 

Otherwise, |xi| [l—| (J^j(Jc) J] < 0 and the infimum is at¬ 
tained at —oo. Hence, the dual function is, 


5(c) = 


Re 


c^yl 


I-^mc||oo < 1 
otherwise. 


(19) 


From (18), \xi\ [l — |(J^j(Jc) J cos= 0 at the infi¬ 


mum, which for every x^ yf 0 yields |(J^^c) J cos (j)i = 1, 


i.e., |(c^J'm)J = 1 and (/>, = 0, as both < 1 


m^ Re (c^y) subject to 


QmxM CmxI 


r.H 

'IxM 


1 


^0, 


M-j 

E] ~ 

i=l 


1, J = 0 
0, J = 1, ■ 


(23) 


1 . 


The number of optimization variables of the dual prob¬ 
lem (23) is (M -I- 1)^/2 equal to half the number of el¬ 


ements of the Hermitian matrix in the inequality con¬ 
straint. Thus, a problem with infinitely many unknown 
parameters (15) is solved over a few optimization vari¬ 
ables. 
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E. Support detection through the dual polynomial 


Strong duality assures that by solving the dual prob¬ 
lem (21), or equivalently Eq. (23), we obtain the min¬ 


imum of the primal problem (15). However, the dual 


problem provides an optimal dual vector, c, but not the 
primal solution, x. Since the corresponding dual polyno¬ 
mial, H{z) = has the properties in Eq. (20), the 

support T of the primal solution x can be estimated by 
locating the angular directions U where the amplitude of 
the dual polynomial is one (i.e., the angular directions at 
the maxima of the beamformed dual vector), 

s i, Vt e i ^ I |7J {z{t))\ < 1, t e T\r. 


Following Sec. m this is done by locating the roots of 
the nonnegative polynomial which lie on the unit circle 
\z\ = 1 (see also Sec. VIII.B), 


M-l 


P{z) = 1 - R{z) = 1 - ^ 

m—— (M—1) 


(25) 


where R{z) = H{z)H{z)^ — \H{z)\‘^ with coefficients 
Cm = m > 0 and r_m = r);^, i.e., the 

autocorrelation of c. 

Note that the polynomial of degree 2(M — 1), 
P+iz)=z^-^P{z) 

1) (26) 


M-l 


= (1^0)^"^-!- ^ rmZ' 

m—— (M —1), m^O 

which has only positive powers of the variable z, has the 
same roots as P{z), besides the trivial root z = 0. Thus, 
the support T of x, i.e., the DOAs of the sources, is 
recovered by locating the roots of P+(z) on the unit circle 
(see Fig.[^, 

f=iu = ^argz, I P+{zi) = 0, \zi\ = l| . (27) 





FIG. 2. (Color online) Support detection through the dual 
polynomial. A ULA is used with M = 21 sensors and 
d/X = 1/2 to localize three sources with support set T = 
[—0.126,0.275,0.67]. (a) The dual polynomial \H{z)\. (b) 
The nonnegative polynomial P{z). (c) The support T is esti¬ 
mated by the angle of the roots, Zi, of P{z) for which \zi\ — 1. 



FIG. 3. (Color online) Grid-free sparse reconstruction. A 
standard ULA is used with M — 21 sensors and d/A = l/2to 
localize three sources (o) at 0 = [—7.2385°, 15.962°, 42.0671°] 
with amplitudes ]x| = [1,0.01,0.6]. (a) The dual polynomial, 
(b) Reconstruction with grid-free CS (*) and CBF. 


F. Reconstruction of the primal solution x 


Once the support is recovered by locating the roots of 
the polynomial in Eq. (26) that lie on the unit circle (27), 
the source amplitudes (the complex weights in Eq. (10)) 
are recovered from, 


XCSd 


= A+y, 


(28) 


where denotes the pseudoinverse of At with columns 
a(U) = eJ2^0/A)[o,-.M-i]^i, fQj. g f 

Figure shows the DOA estimation with grid-free CS 
following the procedure described in this section (see 
App.^for a Matlab implementation). The dual polyno¬ 
mial attains unit amplitude, \H (z)| = 1, at the support 
of the solution, i.e., the DOAs of the existing sources; 
see Fig. I^a). Figure [^b) compares the grid-free CS (28) 


and the CBF (§ reconstruction in DOA estimation. The 
grid-free CS offers very accurate localization, while CBF 
is characterized by low resolution. Moreover, CBF fails 
to detect the weak source at 15.962° since it is totally 
masked by the sidelobes. 


V. MAXIMUM RESOLVABLE DOAS 


The maximum number of resolvable DOAs with grid- 
free CS is determined by the maximum number of roots 
of P+{z) in (26) which can be on the unit circle, jzj = 1. 
Since the coefficients of the polynomial P+{z) are conju¬ 
gate symmetric around the term the roots appear 

in pairs at the same angular direction t;, one inside the 
unit circle, Zin = at radius n < 1, and the 

other outside of the unit circle Zout = = 

l/(zj„)^. This implies that the roots on the unit circle 
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have double multiplicity. The polynomial P+i^) has in 
total 2{M— 1) roots, as determined by its degree. Hence, 
there are at most M —1 (double) roots on the unit circle. 

The necessary condition for the dual polynomial (24) 
to satisfy the condition \H{z)\ < 1 for some t G T, 
thus avoid the non-informative case of a constant dual 
polyn omial , is that the number of sources should not 

exceedi^l^Sl 


K — 

J-^rnn.^r. — 


M - 1 


(29) 


where [-J is the largest integer not greater than the argu¬ 
ment. In other words, at least half of the (paired) M —1 
roots should lie off the unit circle alternating with the 
roots on the unit circle leading to the bound ([2^. 


For positive source amplitudes, Xi G 




the condi¬ 


tion (29) is sufficient and no separation condition is re¬ 
quired for the resolvable sources^^. However, for complex 
amplitudes, Xi G C, the sources are resolved uniquely 
only if t he corresponding DOAs are separated by at 

leastP^I^, 


min \ti — T I 

q.tjGT' 


A 

Md’ 


(30) 



where \ti — tj\ is a wrap-around distance meaning that 
we identify the points —1, 1, in T = [—1,1]. 

The minimum separation condition ( |30[ ) is a conse¬ 
quence of the coherence of the sensing process which is 
related to the beampattern; see Sec. IV.D. in Ref[51 To 
guarantee a well-posed sparse signal reconstruction, it is 
required that the columns of the inverse Fourier opera¬ 
tor J-M, the steering vectors ([^, are sufficiently uncorre¬ 
lated. The continuous formulation (12) implies that ad¬ 


jacent steering vectors are in arbitrarily close directions, 
hence fully coherent. However, the requirement (30) in¬ 


hibits closely spaced (i.e., highly correlated) steering vec¬ 
tors, hence prevents the sparse reconstruction problem 
from being too ill-posed due to coherence. 

Figure shows the reconstruction for the maximum 
number of sources possible. For positive source ampli¬ 
tudes, Xi G M+, the bound (291 suffices to ensure a unique 


solution. Grid-free CS achieves super-resolution even for 
DOAs in general position; see Figs. |^a)-(b). Inserting 
an additional source at 71.81°, thus exceeding the maxi¬ 
mum number of resolvable sources ( |29[ ) , results in a non- 
informative dual polynomial, \H {z) |« 1, for all t G T, 
Fig. s c), and inaccurate reconstruction where only 7 out 
of the 11 sources are resolved. Fig. |^d). For complex 
source amplitudes, Xi G C, an additional constraint (30) 


on the minimum separation of DOAs is required along 
with the bound on the number of sources (29) to ensure 
a unique solution, Figs.|^e)-(f). Violating the minimum 
separation condition, the CS DOA estimation becomes 
extremely ill-posed due to the coherence of the under¬ 
lying steering vectors resulting in inaccurate reconstruc¬ 
tion characterized by the presence of spurious sources, 
Figs.|4];g)-(h). 


FIG. 4. (Color online) Grid-free sparse reconstruc¬ 
tion. A ULA is used with M = 21 sensors and 

d/\ = 1/2 to localize the possible maximum number 

of sources (o), [(M — 1)/2J = 10. (a) The dual 

polynomial and (b) reconstruction with grid-free CS (*) 
and CBF for sources with positive amplitudes, * 10 ,u = 
[0.8, 0.6, 0.9,0.5,1,0.9, 0.1,1,0.4,0.7]. (c) The dual polyno¬ 

mial and (d) reconstruction for 11 sources with positive am¬ 
plitudes, a:ii,R = [xio,K,0.1]. (e) The dual polynomial and (f) 
reconstruction for sources with complex amplitudes, * 10 ,c = 
a:io.R -bi[-1.6, 0.5, -1.3, -2.6, 0.4, -1.2, -1.2, -0.6, -0.5, 0.6], 
separated by the condition (30l. (g) The dual polynomial 


and (h) reconstruction for sources with complex amplitudes, 
a^io.c, but locations violating the condition (30l. 


VI. NON-UNIFORM ARRAYS 

The method is also applicable to non-uniform arrays, 
constructed by randomly choosing sensors from a stan¬ 
dard ULA configuration, by adding an additional con¬ 
straint in the optimization problem (123P3. The ad¬ 


ditional constraint ensures that coefficients of the dual 
polynomial corresponding to inactive sensors on the 
ULA, Cm„uni annihilated. 

The dual problem in a semidefinite programming for¬ 
mulation (231 is augmented with an additional constraint 


and takes the form. 


max Re (c^y) subject to 

c,Q 

M-j , 

Vo ^ ® 


QmxM CmxI 
1 

^IxM ^ 


>1 0, 


(31) 


, M - 1 ’ 
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FIG. 5. (Color online) Grid-free sparse reconstruction, (a) 
A random array constructed by randomly selecting M = 13 
sensors out of a standard ULA with 21 sensors and d/A = 1/2. 
The sources (o) are at = [-32.8881°, 25.2773°, 69.3903°] 
with amplitudes jxj = [0.67, 0.33,1]. (b) The dual polynomial, 
(c) Reconstruction with grid-free CS (*) and CBF. 


Figure shows the DOA estimation with grid-free CS 
and compares it with the CBF reconstruction in the case 
of a random array. Even though CBF performance de¬ 
grades significantly due to the increased sidelobe levels 
introduced by the random array and the strong source 
towards endfire, CS still offers exact reconstruction. 



FIG. 6. (Color online) Grid-free sparse reconstruction. A 
ULA is used with M = 21 sensors and d/\ = 1/2 to localize 
three sources (o) at d = [—19.6942°, 28.3594°, 73.9457°) with 
amplitudes jaij = [0.6, 0.3, 0.3). (a) The dual polynomial, (b) 
Reconstruction with grid-free CS (*) and CBF. The SNR is 
20 dB. 


the support is recovered the source amplitudes are esti¬ 
mated by solving a discrete overdetermined problem (281. 

Figure shows the DOA estimation for three 
sources with grid-free CS when the array measurements 
are contaminated with additive noise (13) such that 
SNR=20 dB. Grid-free CS improves significantly the res¬ 
olution in the reconstruction compared to CBF, even 
though some weak spurious sources appear as artifacts 
due to the noise in the measurements. 


VII. GRID-FREE RECONSTRUCTION WITH NOISE 


The problem of grid-free DOA estimation with CS 
extends to noisy measurements making the framework 
useful for practical applications. Assuming that the 
measurements (131 are contaminated with additive noise 
n G C-^, such that ||n|| 2 < e, the atomic norm minimiza¬ 
tion problem (151 is reformulated a^^, 


min||xm subject to ||y — Fmx \\2 < e. (32) 


To solve the infinite dimensional primal problem ( |32[ ) we 
formulate the equivalent dual problem (see Appendix [ d|, 

max Re (c^y) - e||c ||2 subject to ||-FmcHoo < 1, (33) 

C ^ ^ 


VIII. DOA ESTIMATION WITH POLYNOMIAL 
ROOTING 

Polynomial rooting can increase performance and 
achieve super-resolution in several DOA estimation 
methods, such as the minimum variance distortion¬ 
less response (MVDR) beamformer, the multiple signal 
classification (MUSIC) method and the minimum-norm 
method. All these methods involve the estimation or the 
eigendecomposition of the cross-spectral matrix both in 
their spectral and root version. 

The cross-spectral matrix estimated from L snapshots 
(i.e., observations of y at a particular frequency) is de¬ 
fined as. 


and we replace the infinite-dimensional constraints with 
finite matrix inequalities. 


max 

c.Q 


Re(c^y) —e|jc ||2 sub. to 


QmxM Cjvfxl 
1 

L^ixm 


^0, 


M-j 

Eq 


1, j = 0 

0, J = l,--- ,M-1. 


(34) 


The problem ( |34[ ) is a convex optimization prob¬ 
lem which can be solved efficiently with semidefinite 
programming^ to obtain an estimate for the coefficients, 
c G , of the dual polynomial. The support of the so¬ 
lution, i.e., the DOAs of the existing sources is found by 
locating the points where the dual polynomial has unit 
amplitude following the methodology in Sec jlV.Ej Once 





(35) 


The eigendecomposition of the cross-spectral matrix 
separates the signal and the noise subspaces, 

Cy = U,A«Uf + U„A„U^, (36) 

where Ug comprises the signal eigenvectors, which corre¬ 
spond to the largest eigenvalues Ag, and comprises 
the noise eigenvectors. The signal eigenvectors are in 
the same subspace as the steering vectors (^, while the 
noise eigenvectors are orthogonal to the subspace of the 
steering vectors, thus a(0)^U„ = 0 . 
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A. Spectral version of DOA estimation methods 


MVDRP^ aims to minimize the output power of the 
beamformer under the constraint that the signal from the 
look direction remains undistorted. The MVDR beam- 
former power spectrum is, 


T’mvdr(^) 


1 

"aTeWCyMS)' 


(37) 


MUSICP^ uses the orthogonality between the signal 
and the noise subspace to locate the maxima in the spec¬ 
trum, 


7\iusic(^*) = 


a(0)^U„U«a(0) 


(38) 


The minimum-norm is also an eigendecomposition 
based method but, unlike MUSIC which utilizes all noise 
eigenvectors, it uses a single vector, v = [uq, • • • , vm-i]'^, 
which resides in the noise subspace (compare with the 
dual vector c (20) which resides in the signal subspace) 
such that, 


a(0,)^v = O, i = (39) 


where K is the number of sources. 

All the noise subspace eigenvectors, i.e., the columns 
of U„ have the property in Eq. (39). However, if the 
vector V is chosen as a linear combination of the noise 
subspace eigenvectors the algorithm tends to be more 
robusp™. 


The minimum-norm method selects a vector, v, in the 
noise subspace with minimum ^ 2 -norm and unit hrst ele¬ 
ment, uo = 1. The vector v can be constructed from the 
noise eigenvectors aPl, 


v = U„d^/||d||2, (40) 


where the vector d is the first row of U„. Equivalently, 
the vector v can be constructed from the signal eigenvec¬ 
tors as. 




= U. 




where the vector b is the first row of Ug. 
The minimum-norm spectrum is. 


Hr 


min-norm 


( 0 ) = 


1 


a(0)^vv^a(0) 


(41) 


(42) 


B. Root version of DOA estimation methods 


The root version of the DOA estimation methods is 
based on the fact that for ULAs the null spectrum has 
the form of the trigonometric polynomial in Eq. (B2) with 
uj = 2Tr{d/X) sin0 (since sin0 G [—1,1], then for a stan¬ 
dard ULA w € [—TT, tt]). Thus, evaluating the spectrum 
is equivalent to evaluating the roots of the polynomial on 
the unit circlPI. 


More analytically, let JV(0) = a(0)^^a(0) be the null 
spectrum, such that the spectrum is S(0) = N(0)~^. 
For MVDR, ^ = Cy^ (RefUH p.ll47), for MUSIC, 
= U„U(^ f Ref ITS! p.ll59) and for the minimum-norm 
method, = vv^ (RefUU p.ll63). Then, 

M-l M-1 

sin sine 

m—0 n—0 
M-l 

= ^ (43) 

M-l 

N(z) = ^ 


where ipi = ^mn is the sum of the elements 

of the Hermitian matrix ^ along the lih diagonal and 

2 — gi27r((i/A) sin 0 

The set of DOAs, T, is estimated from the roots 
of the polynomial N{z), or equivalently the polyno¬ 
mial N+{z) = z^~^N{z), which lie on the unit circle, 
as, 

f = |sin6»i = ^ argzj | N+{zi) = 0, \z,\ = l| . (44) 


After the support is recovered, the amplitudes can be 
estimated through an over determined problem as in 

Eq. ([2§. 

Even though the root forms of DOA estimation meth¬ 
ods have, often, more robust performance than the cor¬ 
responding spectral formPl, they require a regular array 
geometry to form a trigonometric polynomial and de¬ 
tect its roots behavior. To achieve a robust estimate of 
the cross-spectral matrix many snapshots are required, 
L > M, i.e., stationary sources. Furthermore, eigende¬ 
composition based methods fail to discern coh erent ar¬ 
rivals. Forward/backward smoothing techniqueP*^ can 
be employed to mitigate this problem and make eigende¬ 
composition based methods suitable for identification of 
coherent sources as well, but they still require a regular 
array geometry and an increased number of sensors. 


IX. EXPERIMENTAL RESULTS 

The high-resolution capabilities of sparse signal recon¬ 
struction methods, i.e., CS for DOA estimation, and the 
robustness of grid-free sparse reconstruction even under 
noisy conditions and with random array configurations 
are demonstrated on ocean acoustic measurements. The 
interest is on single-snapshot reconstruction for source 
tracking and the results are compared with CBF. 

The data set is from the long range acoustic communi¬ 
cations (LRAC) experimenlP^ recorded from 10:00-10:30 
UTC on 16 September 2010 in the NE Pacific and is the 
same as in Ref|5]to allow comparison of the results. The 
data are from a horizontal uniform linear array towed at 
3.5 knots at 200 m depth. The array has M = 64 sen¬ 
sors, with intersensor spacing d = 3 m. The data were 
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acquired with a sampling frequency of 2000 Hz and the 
record is divided in 4 s non-overlapping snapshots. Each 
snapshot is Fourier transformed with 2^^ samples. 

The data are post-processed with CBF and CS on a 
discrete DOA grid [—90°:1°:90°] as well as grid-free CS 
at frequency / = 125 Hz {d/X = 1/4). To facilitate the 
comparison of the results, the grid-free CS reconstruction 
is also presented on the grid [—90°:1°:90°] by rounding 
the estimated DOAs to the closest integer angle and us¬ 
ing the maximum power within each bin. The results are 
depicted in Fig. both with all M — 64 sensors active, 
Figs, [^a)-(d) and by retaining only M = 16 sensors ac¬ 
tive in a non-uniform configuration. Figs, [^e)-(h). Both 
array configurations. Figs. [^a) and [^e), have the same 
aperture thus the same resolution. 

The CBF map in Fig. I^b) indicates the presence 
of three stationary sources at around 45°, 30° and —65°. 
The two arrivals at 45° and 30° are attributed to distant 
transiting ships, even though a record of ships in the area 
was not kept. The broad arrival at —65° is from the tow- 
ship R/V Melville. The CBF map suffers from low reso¬ 
lution and artifacts due to sidelobes and noise. The CS 
reconstruction (§ (e=3.5, Fig. [7]; c)) results in improved 
resolution in the localization of the three sources by pro¬ 
moting sparsity and significant reduction of artifacts in 
the map. The grid-free CS solution (28), Fig. [^d), pro¬ 
vides high resolution and further artifact reduction due 
to polynomial rooting. 

Retaining only 1 /4 of the sensors on the array in a non- 
uniform configuration degrades the resolution of CBF 
due to increased sidelobe levels. Fig. [^f). However, both 
CS on a discrete DOA grid, Figj^g)) and grid-free CS, 
Figj^h), provide high-resolution DOA estimation with¬ 
out a significant reconstruction degradation. 

The single-snapshot processing. Fig. indicates that 
the sources are adequately stationary. Therefore, the 
200 snapshots can be combined to estimate the cross- 
spectral matrix (35) and employ cross-spectral methods 
for DOA estimation. Figure [^a) compares the power 
spectra of MVDR (37), MUSIC ( |M| and the minimum- 
norm method (42) and Fig. [^b) the corresponding root 
versions. 

The root versions of cross-spectral methods, especially 
the root MUSIC and the root minimum-norm method, 
provide improved resolution compared to the correspond¬ 
ing spectral forms. However, the root cross-spectral 
methods require both many snapshots (i.e., stationary 
sources) for a robust estimate of the cross-spectral ma¬ 
trix and uniform arrays. Grid-free CS does not have these 
limitations. 


X. CONCLUSION 

DOA estimation with sensor arrays is a sparse signal 
reconstruction problem which can be solved with com¬ 
pressive sensing (CS). Discretization of the problem in¬ 
volves a compromise between the quality of reconstruc¬ 
tion and the computational complexity, especially for 
high-dimensional problems. Grid-free CS assures that 
the sparsity promoting optimization problem in CS can 
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FIG. 7. (Color online) Data from LRAC. (a) Uniform array 
with M = 64 sensors and the corresponding (b) CBF, (c) 
CS on a discrete grid, [—90°:1°:90°], and (d) grid-free CS 
reconstruction, (e) Non-uniform array with M — 16 sensors 
and the corresponding (f) CBF, (g) CS on a discrete grid and 
(h) grid-free CS reconstruction. 


be solved in the dual domain with semidefinite program¬ 
ming even when the unknowns are infinitely many. Grid- 
free GS achieves high-resolution DOA estimation through 
the polynomial rooting method. 

In contrast to established DOA estimation meth¬ 
ods, CS provides high-resolution acoustic imaging even 
with non-uniform array configurations and robust per¬ 
formance under noisy measurements and single-snapshot 
data. Finally, the grid-free CS has the same performance 
both with coherent and incoherent, stationary or moving 
sources while other DOA estimation methods based on 
polynomial rooting fail to discern coherent arrivals and 
have degraded resolution for moving sources as they re¬ 
quire many snapshots. 
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FIG. 8. (Color online) Data from LRAC, combining the 200 
snapshots to estimate the cross-spectral matrix and process¬ 
ing with MVDR, MUSIC and the minimum-norm method, 
(a) Spectral version and (b) root version. The ULA with 
M = 64 sensors and d/A = l/4is used. 


APPENDIX A: CONVEX OPTIMIZATION PROBLEMS 

This section summarizes the basic notions and formu¬ 
lations encountered in convex optimization problems, as 
presented analytically in ReflMl 


1. Primal problem 

A generic optimization problem has the form, 
min /o(x) 

X 

subject to /i(x) <0, i = 1, • • • , TO (Al) 

/ij(x) = 0, j = 1, - • • 

where x S is the optimization variable, the function 
/o : —>• M is the objective (or cost) function, the func¬ 
tions fi : —>■ K are the inequality constraint functions 

and the functions hj : —>■ C are the equality con¬ 

straint functions. The optimization problem is con¬ 
vex when /o, • ■ ■ > /m are convex functions and hi, • • • ,hq 
are affine (linear) functions. 

The set of points for which the objective and all con¬ 
straint functions in Eq. (A1) are defined is called the 


domain of the optimization problem, 

m q 


I? = Pi dom/i n Pi domhj. 

1=1 


(A2) 


z=0 


A point X G I? is called feasible if it satisfies the con¬ 
straints in Eq. (All. 


The optimal value p* of the optimization prob¬ 
lem (Al I, achieved at the optimal variable x*, is. 


(A3) 


p* = inf {/o(x) I /i(x) < 0, hj{x) = 0} 
= {/o(x*) I /*(x*) < 0, hj{x*) = 0} , 
for alH = 1, • • • , TO and j = 1, • • • ,q. 

2. The Lagrangian 


The Lagrangian, L, of an optimization problem is 
obtained by augmenting the objective function with a 


weighted sum of the constraint functions. The La¬ 
grangian of the generic optimization problem (Al) is , 


L(x, A, zz) =/o(x)-f ^ Ai/*(x)-fRe 


Z=1 


q 

E 

1=1 


Vihi{x) 


(A4) 


where Xi is the Lagrange multiplier associated with the 
zth inequality constraint, /i(x) < 0, and Vj is the La¬ 
grange multiplier associated with the jth equality con¬ 
straint, hj{x.) = 0. The vectors A G K"* and u G are 
the dual variables of the problem (Al). 


3. The dual function 


The dual function of the problem (Al) is the minimum 
value of the Lagrangian ( A4 ) over x G I? for A G K"* and 
jz G C«, 


g(A,iz) = inf L(x, A,zz). 


(AS) 


Since the dual function is the pointwise infinum of a fam¬ 
ily of affine functions of (A, iz), it is concave, even when 
the problem (|A1[) is not convex. 


The dual function (AS) yields lower bounds on the op¬ 
timal value p* 


for any A ^ 0 (where ^ represents 
componentwise inequality) and any iz. 


since g{X, zz) = inf L(x, A 


5(A,iz)<p*, (A6) 

zz) < L(x,A,zz) < /o(x) for 


every feasible point x. 


4. Dual problem 


The dual function (AS) gives a lower bound on the op¬ 


timal value p* of the optimization problem (All, which 


depe nds on the dual variables (A, zz) with A ^ 0; see 
Eq. (A6|. The best lower bound, i.e. the lower bound 


with the greatest value, is obtained through the optimiza¬ 
tion problem. 


max g{\, zz) subject to A ^ 0, 

A,t> 


(A7) 


which is the dual problem to the optimization prob¬ 
lem (Al I. 


The dual problem (A7) is a convex optimization prob¬ 


lem, since the objective function to be maximized is 
concave and the constraints are convex, irrespectively 
whether the primal problem (Al) is convex or not. 


5. Weak duality 


The optimal value d* of the dual problem (A71, 
achieved at the dual optimal variables (A*,zz*) is. 


d* = sup{ 5 (A,zz) I A ^ 0} 
= {g(A*,zz*) I A* ^0}. 


(AS) 
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The dual maximum d* is the best lower bound on the 


minimum of the primal problem (A3), that can be ob¬ 


tained from the Lagrange dual function. The inequality, 
d* <p*, (A9) 


holds even if the primal problem ( |Al| is non-convex and 
is called weak duality. 

The non-negative difference p* — d* is called the duality 
gap for the optimization problem (A1), since it gives the 


gap between the minimum of the primal problem and the 
maximum of the dual problem. 


6. Slater’s condition and strong duality 

When the duality gap, p* — d*, is zero, strong duality 
holds characterized by the equality. 


APPENDIX B: BOUNDED TRIGONOMETRIC 
POLYNOMIALS 

This section presents useful results for bounded 
trigonometric polynomials and their roots as presented 
in RefsjSZl 


1. Trigonometric polynomials 

Let a(w) = he a L x 1 basis 

vector for trigonometric polynomials of degree L — 1 with 
u! G [—TT, tt]. a (causal) trigonometric polynomial can be 
written in terms of the basis vector as, 

L-l 

H{uj) = = a(w)^h, (Bl) 

1=0 


d* =p*. 


(AlO) 


Strong duality holds when the optimization prob¬ 
lem (A1) is convex and there exists a strictly feasible 
point, i.e., the inequality constraints hold with strict 
inequalities. The constraint qualification which implies 
strong duality for convex problems is called Slater’s con¬ 
dition. 


/i(x) <0, i = !,■■■ ,m, 
Aqxjvx = y. 


(All) 


When the primal problem is convex and Slater’s condi¬ 
tion holds there exist a dual feasible (A*,zz*) such that 
g{\*, u*) = d* = p*, i.e., the optimal value of the primal 
problem can be obtained by solving the dual problem. 

The Slater’s condition holds also with a weaker con¬ 
straint qualification, when some of the inequality con¬ 
straint functions, /i,-- - ,/fe, are affine (instead of con¬ 
vex). 


/*(x) < 0, i = I,- - • ,fc, 

/j(x) < 0, i = fc -I- 1, • • • , m, (A12) 

Agxwx = y. 


The weaker constraint qualifications (A12) imply that 
strong duality reduces to feasibility when both the in¬ 
equality and the equality constraints are linear. 


7. Schur complement 


where h = [Hq,--- G is the vector of the 

polynomial coefficients. 


2. Nonnegative trigonometric polynomials 


Let R{oj) = \H(uj)\'^ = H{uj)H{uj)^. From (Bl|, 
the nonnegative trigonometric polynomial R{uj) has the 
form. 


L-l 

R{uj) = Y (B2) 

fe=-(L-l) 

L-l-k 

where ^ for fc > 0 and r_fe = r^, i.e., the 

1=0 

coefficients are conjugate symmetric thus R{u}) is Hermi- 
tian. Equivalently, the coefficients can be calculated 
as the sum of the diagonal elements of the autocor¬ 
relation matrix QixL = hh^ as, 


L-k 

i=l 


3. Bounded trigonometric polynomials 

Let two polynomials H{uj) and fulfill the inequal¬ 
ity, 


\H{oj)\ < \B{uj)\, "iw G [-7r,7r], (B4) 


Let X be a square Hermitian matrix partitioned as, 


A B 
C ’ 


(A13) 


where A is also square Hermitian. If det A 0 then the 
matrix, 

S = C-B^A-iB, (A14) 

is called the Schur complement of A in X. 

A useful property related to the Schur complement is 
that if A 0 then X ^ 0 if and only if S ^ 0. 


which implies |iL(a;)p < |i?(w)p, Vw G [—7r,7 r]. D efining 
i?_fy(a;) = \H{uj)\^ and Rb{<jj) = \B{uj)\^ as in (B2), yields 
Rnii^) < RBioj)- From Lemma 4.23 iiP^, it!_fy(w) < 
Rb{^) implies Q/y ^ Qs, where Q/y = hh^ and Qs = 
bb^ are the autocorrelation matrices of the coefficient 
vectors h = [/iq, • • • , and b = [6o, • • • , of 

the polynomials H{uj) and Bjuj) respectively. Through a 
Schur complement (see Sec. A.7), Qb — hl“^h^ ^ 0 is 
equivalent to semidefinite matrix , 


Qb 

y^H 

L*^1xL 


hixi 

1 


^ 0 . 


(B5) 
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Let the polynomial Hiui) have amplitude uniformly 
bounded for all uj G [—such that, \H{uj)\ < 7 , where 
7 € K+ is a given positive real number. As a special case 
of the results for bounded trigonometric polynomials in 
Eqs. (B4), (B5), with |i?(a;)| = 7 , Theorem 4.24 and 
corollary 4.25 iipZl states that the inequality |i^(a;)| < 7 
can be approximated by two linear matrix inequalities, 


QlxL h^xl 

hfxL 1 


L-J 


10 , 


E ^ ^ 


^ j = 0 

0, J = 1, 


(B6) 


,L-1. 


The latter constraint follows from the autocorrelation 
matrix of the constant polynomial Rb{oj) = 7 ^- 

The results for bounded trigonometric polynomials can 
be used in relation to the £oo-norm, since setting an up¬ 
per bound for the maximum amplitude of a polynomial 
implies that the polynomial has amplitude uniformly 
bounded for all w G [—7r,7r], 

|jiJ||oo = max |iL(a;)| < 7 , 

Uie[—!7,TV] 

\H{uj)\ <^,yujG [—7r,7r]. 


4. Roots of real nonnegative trigonometric polynomials 

For a bounded trigonometric polynomial \H(u})\ < 1, 
we can construct a polynomial, 

P{uj) = l-\H{uj)\‘^ = 1-R{u}), (B8) 


TABLE I. Matlab code for Sec. US 


Given y G , d, A 



Solve dual problem with 

CVXP, Eq. 

(231 

1 

cvx_solver sdpt3 



2 

cvx_begin sdp 



3 

variable S{M + 1, M 

-b 1) hermitian 

4 

s>=o-, 



5 

S{M + 1,M + 1) == 

1; 


6 

trace(5') == 2; 



7 

for j = 1 : M — 1 



8 

sum(diag(S', j)) == 


j,M-bl); 

9 

end 




10 

11 

12 


maximize(real(S'(l : M, M + 1)' * y)) 
cvx_end 

c = ^(l : + 


Find the roots of P+, Eq |26l 


13 

14 

15 


r = conv(c,fiipud(conj(c))); 
r{M) = l-r(M); 
roots P = roots(r); 


Isolate roots on the unit circle, Eq. (271 


16 

17 

18 
19 


roots_uc = roots_P(abs(l-abs(roots_P))< le — 2); 
[auxjind] =sort (real(roots_uc)); 
roots_uc = roots_uc(ind); 

t = angle(roots uc(l:2:end))/(2 * pi * d/lam6da); 


Amplitude estimation, Eq. (281 


20: A_T = exp(li * 2* pi * d/lambda*[0:{M — 1)]' * t'); 
21: x CS dual = A T\i/; 


which is by definition real-valued and nonnegative, thus 
it cannot have single roots on the unit circle. The de¬ 
gree of the polynomial P{uj) is 2{L — 1). Therefore, the 
polynomial P(w) has at most L — 1 distinct roots on the 
unit circle. At a root, luq, we have P(ojo) = 0 and subse¬ 
quently |iJ(a;o)| = 1. 

APPENDIX C: IMPLEMENTATION IN MATLAB 

The algorithm in Tabl e |l| fo r the implementation of the 
method described in Sec.|l^is an adaptation of the code 
by Fernandez-Granda in ReflTTl 


APPENDIX D: DUAL PROBLEM WITH NOISE 


constraints, 

L(x,c,^) = 

^ ' fD2') 

||a:|U + Re [c^ (y - Pmx - n)] -|- ^ (n^n - e^) , 

where c G are the dual variables related to the equal¬ 
ity constraints, y — Pmx — n = 0, and ^ G K+ is a 
Lagrange multiplier related to the inequality constraint, 
|ln|l2-e<0. 

The dual function g{c,^) is the infimum of the La- 
grangian, L(x,c,^), over the optimization variable x, 

5(c,0 = inf L{x,c,C) 

X 

= Re [c^y — c^n] -f ^ (n^n — e^) -b (D3) 
+ inf (||a;|U - Re [c^P'Ma^]) ■ 


In the case that the measurements are contami¬ 
nated with additive noise n G such that ||n ||2 < e. 


the primal problem of atomic norm minimization (151 is 
reformulated to the problem (32) or equivalently. 


minllccllxt subject to 


y = Pmx + n, 
l|n ||2 < e. 


(Dl) 


The Lagrangian for (Dl) is formulated by augment¬ 


ing the objective function with a weighted sum of the 


Minimizing over the unknown noise n G 

dg{c,0 


On 


= —c -b 2^n = 0, 


(D4) 


yields the optimal noise vector, Hq = c/(2^). The dual 
function evaluated at Oq is, 

-b inf (||x|U - Re [c^Pmx] ) . 
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(D6) 


Further, maximizing over the dual variable 

^ _ 2 ^ n 

4^2 e 

we obtain the optimal value for the dual variable = 
Ilc|l2/(2e). 

Finally, the dual function evaluated at the optimal val¬ 
ues Ho and becomes. 


5(c)|no.5o = Re [c^y] - e||c|| 2 -f 

-finf (||a;|U - Re [c^J'Ma;]) 


(D7) 


and the dual problem is formulated by maximizing the 
dual function, g(c)|no,4oi over the dual variables c G 
similarly to the process detailed in Sec. |IV.C[ 

max5(c)|„^_5^ = 


max Re [c^y] - e||c|j 2 subject to ||-F^c||oo < 1- 


(D8) 
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